Analysis

Виконавці: Данило Кравчук (завантаження даних, DESeq, Функц. аналіз) Лев Кравчук (PCA, Функц. аналіз)

GitHub

Мета

Для аналізу диференційної експресії генів було обрано клітини лінії А549 як біологічний об’єкт. Клітини лінії A549 — це клітини недрібноклітинної (NSCLC) аденокарциноми легень.

Завантаження даних

Дані було отримано з бази даних GEO (Gene Expression Omnibus). Завантажені дані описують експеримент з секвенуванням РНК (RNA-sequencing) на опромінених рентгенівським випромінюванням і контрольних клітинах A549.

Клітини A549 були піддані одноразовій дозі рентгенівського опромінення 2Гр або 4Гр за допомогою приладу X-RAD 160-225. Після опромінення клітини A549 (з одноразовою дозою опромінення), а також неопромінені клітини A549 були інкубовані ще 48 годин. Після цього було зібрано загальну РНК для високопродуктивного секвенування. Загальний дизайн експерименту: дослідження 3 умов по 3 репліки. Три групи: група 0Гр (контроль), група 2Гр та група 4Гр.

Відповідна стаття: doi: 10.1038/s41419-022-04561-x, PMID: 35190532.

Дані експерименту та метадані

Метадані зразків

kable(rse_colData_f |> 
  select(external_id, sra.study_title, sra.sample_title) |>  head(), format = "html", table.attr = "class='table table-striped'")

	external_id	sra.study_title	sra.sample_title
SRR8371688	SRR8371688	RNA-sequencing in irradiated and normal A549 cells.	0GY-1
SRR8371689	SRR8371689	RNA-sequencing in irradiated and normal A549 cells.	0GY-2
SRR8371690	SRR8371690	RNA-sequencing in irradiated and normal A549 cells.	0GY-3
SRR8371691	SRR8371691	RNA-sequencing in irradiated and normal A549 cells.	2GY-1
SRR8371692	SRR8371692	RNA-sequencing in irradiated and normal A549 cells.	2GY-2
SRR8371693	SRR8371693	RNA-sequencing in irradiated and normal A549 cells.	2GY-3

Метадані генів

kable(head(rse_rowData_f), format = "html", table.attr = "class='table table-striped'")

	source	type	bp_length	phase	gene_id	gene_type	gene_name	level	havana_gene	tag
ENSG00000278704.1	ENSEMBL	gene	2237	NA	ENSG00000278704.1	protein_coding	BX004987.1	3	NA	NA
ENSG00000277400.1	ENSEMBL	gene	2179	NA	ENSG00000277400.1	protein_coding	AC145212.2	3	NA	NA
ENSG00000274847.1	ENSEMBL	gene	1599	NA	ENSG00000274847.1	protein_coding	AC145212.1	3	NA	NA
ENSG00000276256.1	ENSEMBL	gene	2195	NA	ENSG00000276256.1	protein_coding	AC011043.1	3	NA	NA
ENSG00000278198.1	ENSEMBL	gene	1468	NA	ENSG00000278198.1	protein_coding	AC011043.2	3	NA	NA
ENSG00000273496.1	ENSEMBL	gene	1468	NA	ENSG00000273496.1	protein_coding	AC011841.1	3	NA	NA

Матриця експресії

kable(rse_assay_f |> 
  arrange(desc(SRR8371688)) |> head(), format = "html", table.attr = "class='table table-striped'")

	SRR8371688	SRR8371689	SRR8371690	SRR8371691	SRR8371692	SRR8371693	SRR8371694	SRR8371695	SRR8371696
ENSG00000198804.2	11438923	12624134	9876942	11286817	11500584	7659250	5559007	6564597	9430680
ENSG00000115414.18	6535629	7642890	5454817	6763706	5559684	4748784	3533072	3375268	4997642
ENSG00000198886.2	6206643	7357551	5660816	5944139	5889170	3652051	2882935	3346084	4630646
ENSG00000165092.12	5324460	5907729	4752635	6099076	6003428	4324184	2738698	2820161	4443411
ENSG00000198712.1	3950989	4529637	3600388	4185883	4218314	2704965	2025084	2605506	3313521
ENSG00000198786.2	3844190	4525495	3422576	3768348	3622434	2426284	1887319	2046555	3086360

Попередня візуалізація даних

Побудовано графік щільності для окремого зразка (SRR8371688) та гістограму частоти значень експресії

# Створення таблиці метаданих
colData_new <- data.frame(
  row.names = colnames(rse_filtered),
  condition = factor(c(rep("0Gy", 3), rep("2Gy", 3), rep("4Gy", 3)))
)

# Матриця counts
countData <- assay(rse_filtered)

# Побудова графіку щільності
p1 <- ggplot(countData, aes(x = SRR8371688 )) +
  geom_density(color = "darkblue", fill = "lightblue") +
  xlim(0, 10e4)

count_summary <- countData |>
  as.vector() |>
  data.frame(counts = _) |>
  group_by(counts) |>
  summarise(frequency = n(), .groups = 'drop') |>
  arrange(counts)


# Створення гістограми частоти значень
p2 <- ggplot(count_summary[1:200,], aes(x = counts, y = frequency, fill = frequency)) +
  geom_bar(stat = "identity", width = 0.7) +
  scale_y_log10(labels = scales::comma) +
  labs(title = "Histogram of Counts vs Frequency (Logarithmic Scale)",
       x = "Counts",
       y = "Frequency (log10)") +
  scale_fill_gradient(low = "#AA05BA", high = "#56B1F7") +
  theme_minimal()
#p1
#p2
counts_info <- p1 + p2

ggsave("./figs/counts_info.png", plot = counts_info, width = 16, height = 6, dpi = 300)
knitr::include_graphics("./figs/counts_info.png")

Аналіз диференційної експресії генів

DESeq

Для аналізу диференційної експресії генів було застосовано DESeq Wald

Коригування p.adj

Додатково було проведено коригування на хибний позитивний результат p.adj (FDR). Нижче наведено відповідний графік

res_wald_2Gy_vs_0Gy <- results(dds_wald, contrast = c("condition", "2Gy", "0Gy"))
res_wald_4Gy_vs_0Gy <- results(dds_wald, contrast = c("condition", "4Gy", "0Gy"))
res_wald_4Gy_vs_2Gy <- results(dds_wald, contrast = c("condition", "4Gy", "2Gy"))

hist_data_before <- bind_rows(
  data.frame(padj = res_wald_2Gy_vs_0Gy$padj, Comparison = "2Gy vs 0Gy"),
  data.frame(padj = res_wald_4Gy_vs_0Gy$padj, Comparison = "4Gy vs 0Gy"),
  data.frame(padj = res_wald_4Gy_vs_2Gy$padj, Comparison = "4Gy vs 2Gy")
)

res_wald_2Gy_vs_0Gy$padj <- p.adjust(res_wald_2Gy_vs_0Gy$pvalue, method = "fdr")
res_wald_4Gy_vs_0Gy$padj <- p.adjust(res_wald_4Gy_vs_0Gy$pvalue, method = "fdr")
res_wald_4Gy_vs_2Gy$padj <- p.adjust(res_wald_4Gy_vs_2Gy$pvalue, method = "fdr")

hist_data_after <- bind_rows(
  data.frame(padj = res_wald_2Gy_vs_0Gy$padj, Comparison = "2Gy vs 0Gy"),
  data.frame(padj = res_wald_4Gy_vs_0Gy$padj, Comparison = "4Gy vs 0Gy"),
  data.frame(padj = res_wald_4Gy_vs_2Gy$padj, Comparison = "4Gy vs 2Gy")
)

deseq_Wald_padj_b <- ggplot(hist_data_before, aes(x = padj, fill = Comparison, color = Comparison)) +
  geom_density(alpha = 0.2, adjust = 1.5) + # Криві густини
  scale_fill_manual(values = c("2Gy vs 0Gy" = "red", "4Gy vs 0Gy" = "blue", "4Gy vs 2Gy" = "forestgreen")) +
  scale_color_manual(values = c("2Gy vs 0Gy" = "red", "4Gy vs 0Gy" = "blue", "4Gy vs 2Gy" = "forestgreen")) +
  labs(title = "Розподіл коригованих p-значень", x = "p-значення", y = "Щільність") +
  geom_vline(xintercept = 0.05, linetype = "dashed", color = "black") +
  annotate("text", x = 0.05, y = Inf, label = "0.05", vjust = 1.5, hjust = -0.1, color = "black") +
  theme_minimal()

deseq_Wald_padj_a <- ggplot(hist_data_after, aes(x = padj, fill = Comparison, color = Comparison)) +
  geom_density(alpha = 0.2, adjust = 1.5) + # Криві густини
  scale_fill_manual(values = c("2Gy vs 0Gy" = "red", "4Gy vs 0Gy" = "blue", "4Gy vs 2Gy" = "forestgreen")) +
  scale_color_manual(values = c("2Gy vs 0Gy" = "red", "4Gy vs 0Gy" = "blue", "4Gy vs 2Gy" = "forestgreen")) +
  labs(title = "Розподіл коригованих p-значень", x = "p-значення", y = "Щільність") +
  geom_vline(xintercept = 0.05, linetype = "dashed", color = "black") +
  annotate("text", x = 0.05, y = Inf, label = "0.05", vjust = 1.5, hjust = -0.1, color = "black") +
  theme_minimal()

deseq_Wald_padj <- deseq_Wald_padj_b + deseq_Wald_padj_a
ggsave("./figs/deseq_Wald_padj.png", plot = deseq_Wald_padj, width = 16, height = 6, dpi = 300)
knitr::include_graphics("./figs/deseq_Wald_padj.png")

Volcano plot

Для візуалізації у зміні експресії генів для випадку опромінення рентгеном 2Гр (порівняння з контролнм) та 4Гр (порівняння з контролнм) побудовано volcano plot. Окремо проводилось порівняння експресії генів при опроміненні 2Гр та 4Гр.

# png("./figs/deseq_Wald.png", width = 2400, height = 800)
# par(mfrow = c(1, 3), mar = c(5, 5, 4, 2))

# plotMA(res_wald_2Gy_vs_0Gy, ylim=c(-30,30), main="MA-plot: 0Gy vs 2Gy")
# plotMA(res_wald_4Gy_vs_0Gy, ylim=c(-30,30), main="MA-plot: 0Gy vs 4Gy")
# plotMA(res_wald_4Gy_vs_2Gy, ylim=c(-30,30), main="MA-plot: 2Gy vs 4Gy")
# dev.off()
# knitr::include_graphics("./figs/deseq_Wald.png")

res_wald_df_2Gy_vs_0Gy <- as.data.frame(res_wald_2Gy_vs_0Gy)
res_wald_df_2Gy_vs_0Gy$significant <- res_wald_df_2Gy_vs_0Gy$padj < 0.05

res_wald_df_4Gy_vs_0Gy <- as.data.frame(res_wald_4Gy_vs_0Gy)
res_wald_df_4Gy_vs_0Gy$significant <- res_wald_df_4Gy_vs_0Gy$padj < 0.05

res_wald_df_4Gy_vs_2Gy <- as.data.frame(res_wald_4Gy_vs_2Gy)
res_wald_df_4Gy_vs_2Gy$significant <- res_wald_df_4Gy_vs_2Gy$padj < 0.05

res_wald_df_2Gy_vs_0Gy$comparison <- "2Gy vs 0Gy"
res_wald_df_4Gy_vs_0Gy$comparison <- "4Gy vs 0Gy"
res_wald_df_4Gy_vs_2Gy$comparison <- "4Gy vs 2Gy"

res_wald_df_2Gy_vs_0Gy <- as.data.frame(res_wald_df_2Gy_vs_0Gy)  |> 
  mutate(gene_id = rownames(res_wald_df_2Gy_vs_0Gy), comparison = "2Gy vs 0Gy")  |> 
  left_join(gene_annotations, by = "gene_id")

# Оновлення res_df_4Gy_vs_0Gy
res_wald_df_4Gy_vs_0Gy <- as.data.frame(res_wald_df_4Gy_vs_0Gy)  %>% 
  mutate(gene_id = rownames(.), comparison = "4Gy vs 0Gy")  |> 
  left_join(gene_annotations, by = "gene_id")

# Оновлення res_df_4Gy_vs_2Gy
res_wald_df_4Gy_vs_2Gy <- as.data.frame(res_wald_df_4Gy_vs_2Gy)  |> 
  mutate(gene_id = rownames(res_wald_df_4Gy_vs_2Gy), comparison = "4Gy vs 2Gy")  |> 
  left_join(gene_annotations, by = "gene_id")

# Побудова volcano plot
plot_2Gy_vs_0Gy_W <- EnhancedVolcano(res_wald_df_2Gy_vs_0Gy,
                                   lab = res_wald_df_2Gy_vs_0Gy$gene_name,
                                   x = 'log2FoldChange',
                                   y = 'padj',
                                   xlim = c(-5, 5),
                                   ylim = c(0, 10),
                                   pCutoff = 0.05,
                                   FCcutoff = 1.5,
                                   pointSize = 5,
                                   labSize = 6,
                                   col = c('grey30', 'forestgreen', 'royalblue', 'red2'),
                                   title = '2Gy vs 0Gy',
                                   legendPosition = 'top')

plot_4Gy_vs_0Gy_W <- EnhancedVolcano(res_wald_df_4Gy_vs_0Gy,
                                   lab = res_wald_df_4Gy_vs_0Gy$gene_name,
                                   x = 'log2FoldChange',
                                   y = 'padj',
                                   xlim = c(-5, 5),
                                   ylim = c(0, 10),
                                   pCutoff = 0.05,
                                   FCcutoff = 1.5,
                                   pointSize = 5,
                                   labSize = 6,
                                   col = c('grey30', 'forestgreen', 'royalblue', 'red2'),
                                   title = '4Gy vs 0Gy',
                                   legendPosition = 'top')

plot_4Gy_vs_2Gy_W <- EnhancedVolcano(res_wald_df_4Gy_vs_2Gy,
                                   lab = res_wald_df_4Gy_vs_2Gy$gene_name,
                                   x = 'log2FoldChange',
                                   y = 'padj',
                                   xlim = c(-5, 5),
                                   ylim = c(0, 10),
                                   pCutoff = 0.05,
                                   FCcutoff = 1.5,
                                   pointSize = 5,
                                   labSize = 6,
                                   col = c('grey30', 'forestgreen', 'royalblue', 'red2'),
                                   title = '4Gy vs 2Gy',
                                   legendPosition = 'top')
volc_wald_plot <- plot_2Gy_vs_0Gy_W + plot_4Gy_vs_0Gy_W + plot_4Gy_vs_2Gy_W + plot_layout(nrow = 1)

#print(volc_wald_plot)
ggsave("./figs/volc_wald.png", plot = volc_wald_plot, width = 20, height = 6, dpi = 300)

2Gy vs 0Gy

ggsave("./figs/volc_plot_2Gy_vs_0Gy.png", plot = plot_2Gy_vs_0Gy_W, width = 20, height = 20, dpi = 300)
knitr::include_graphics("./figs/volc_plot_2Gy_vs_0Gy.png")

4Gy vs 0Gy

ggsave("./figs/volc_plot_4Gy_vs_0Gy.png", plot = plot_4Gy_vs_0Gy_W, width = 20, height = 20, dpi = 300)
knitr::include_graphics("./figs/volc_plot_4Gy_vs_0Gy.png")

4Gy vs 2Gy

ggsave("./figs/volc_plot_4Gy_vs_2Gy.png", plot = plot_4Gy_vs_2Gy_W, width = 20, height = 20, dpi = 300)
knitr::include_graphics("./figs/volc_plot_4Gy_vs_2Gy.png")

Upset

Для візуалізації перекриття між списками генів, які є статистично значущими (з коригованим p-значенням < 0.05) та мають великий зсув логарифмічних змін (log2FoldChange > 1.5 для підвищеної експресії та log2FoldChange < -1.5 для зниженої експресії) у різних порівняннях, використали upset plot.

res_2Gy_0Gy_up <- res_wald_df_2Gy_vs_0Gy |> as.data.frame() |> filter(pvalue < 0.05 & log2FoldChange > 1.5)
res_2Gy_0Gy_down <- res_wald_df_2Gy_vs_0Gy  |> as.data.frame() |> filter(pvalue < 0.05 & log2FoldChange < -1.5)
res_4Gy_0Gy_up <- res_wald_df_4Gy_vs_0Gy  |> as.data.frame() |> filter(pvalue < 0.05 & log2FoldChange > 1.5)
res_4Gy_0Gy_down <- res_wald_df_4Gy_vs_0Gy  |> as.data.frame() |> filter(pvalue < 0.05 & log2FoldChange < -1.5)
res_4Gy_2Gy_up <- res_wald_df_4Gy_vs_2Gy |> as.data.frame() |> filter(pvalue < 0.05 & log2FoldChange > 1.5)
res_4Gy_2Gy_down <- res_wald_df_4Gy_vs_2Gy  |> as.data.frame() |> filter(pvalue < 0.05 & log2FoldChange < -1.5)

venn_list <- list("2Gy vs 0Gy(up)" = res_2Gy_0Gy_up$gene_name,
                  "2Gy vs 0Gy(down)" = res_2Gy_0Gy_down$gene_name,
                  "4Gy vs 0Gy(up)" = res_4Gy_0Gy_up$gene_name,
                  "4Gy vs 0Gy(down)" = res_4Gy_0Gy_down$gene_name,
                  "4Gy vs 2Gy(up)" = res_4Gy_2Gy_up$gene_name,
                  "4Gy vs 2Gy(down)" = res_4Gy_2Gy_down$gene_name)

ven <- venndetail(venn_list)

png("./figs/upset_diagram.png", width = 1600, height = 1600, res = 300)
plot(ven, type = "upset")
dev.off()

Аналіз головних компонент (PCA)

Наступним кроком було проведено PCA аналіз.

Scree plot

Червона лінія відображає кумулятивну суму варіацій, яка показує, що 80% варіації даних припадає на п’ять головних компонент. При цьому PC1 пояснює близько 50% варіації. Для того, щоб визначити кількість головних компонент, які варто зберегти для аналізу були використані метод Горна та elbow, які зображені на графіку штриховою лінією. У більшості випадків знайдені значення не збігатимуться. Це пов’язано з тим, що визначення правильної кількості головних компонент є складним завданням і визначення “правильної” кількості кластерів не існує.

Biplot

Він комбінує інформацію про зразки і змінні в одному графіку. По осях зображені дві головні компоненти: PC1 і PC2, кожна з яких пояснює відсоток варіації в даних: PC1 — 44.98%, PC2 — 12.94%. Точки відповідають різним спостереженням за умов випромінювання X-хвилями 0Гр, 2Гр, 4Гр, які приблизно кластеризуються відповідно умовам впливу. Вектори представляють гени, які мають найбільший вплив на кожну компоненту. Напрям стрілки вказує, як експресія гену впливає на компоненти, а довжина — наскільки сильний цей вплив.

Щоб перевірити чи є “правильна” кластеризація експериментів, був проведений кластерний аналіз методом k-means з 5-ма головними компонентами, адже саме стільки компонент пояснюють більшу варіацію даних.

Щоб представити взаємодію між різними головними компонентами PC1, PC2, PC3, PC4, PC5, побудовано Pairs plot. На всіх графіках видно кластеризацію даних, але найліпше це спостерігається на парах PC1 - PC2-5 та PC2 - PC3-4.

Loadings plot

Він показуює, наскільки кожен ген впливає на кожну головну компоненту. Це значення визначає, наскільки цей ген сприяє поясненню варіації в даних. Гени з високими позитивними або негативними завантаженнями мають найбільший внесок у головні компоненти та позначені фіолетовим та жовтим кольором відповідно.

vst_data <- assay(vst(dds_wald))

pca_result <- PCAtools::pca(vst_data, metadata = colData_new, removeVar = 0.1)

horn <- parallelPCA(vst_data)
elbow <- findElbowPoint(pca_result$variance)
var_explained <- which(cumsum(pca_result$variance)> 80)[1]

pscree <- PCAtools::screeplot(pca_result,
    axisLabSize = 18,
    titleLabSize = 22,
    vline = c(horn$n, elbow, var_explained), hline = 85,
    returnPlot = FALSE) +
    geom_label(aes(x = horn$n, y = 60,
      label = 'Horn\'s', vjust = -1, size = 8)) +
    geom_label(aes(x = elbow, y = 60,
      label = 'Elbow method', vjust = -3, size = 8)) +
    geom_label(aes(x = var_explained , y = 85,
      label = '80% explained variation', vjust = -1, size = 8)) 


ppairs <- PCAtools::pairsplot(pca_result,
    hline = 0, vline = 0,
    pointSize = 2,
    gridlines.major = FALSE, gridlines.minor = FALSE,
    colby = 'condition',
    plotaxes = TRUE,
    margingaps = unit(c(0.01, 0.01, 0.01, 0.01), 'cm'),
    returnPlot = FALSE)


# Створюємо нові назви для рядків, зіставляючи gene_id з gene_name
new_row_names <- gene_annotations$gene_name[match(rownames(pca_result$loadings), gene_annotations$gene_id)]

# Якщо є дублікат імен, додамо індекси для унікальності
new_row_names <- make.unique(new_row_names)

# Перезаписуємо назви рядків у pca_result$loadings
rownames(pca_result$loadings) <- new_row_names


pbiplot <- PCAtools::biplot(pca_result,
    axisLabSize = 18,
    # loadings parameters
    showLoadings = TRUE,
    sizeLoadingsNames = 4,
    colLoadingsNames = 'red4',
    #points
    lab = NULL,
    colby = 'condition',
    colLegendTitle = 'Condition',
    colkey = c('0Gy' = '#228B22', '2Gy' = '#7B68EE', '4Gy' = '#FF8C00'),
    pointSize = 5,
    #plotview
    hline = 0,
    vline = c(-5, 0, 5),
    vlineType = c('dotdash', 'solid', 'dashed'),
    gridlines.major = FALSE, gridlines.minor = FALSE,
    legendPosition = 'right', 
    legendLabSize = 16,
    legendIconSize = 8.0,,
    subtitle = 'PC1 versus PC2',
    caption = '5 PCs ≈ 80%',
    returnPlot = FALSE)

ploadings <- PCAtools::plotloadings(pca_result,
    axisLabSize = 18,
    subtitle = 'PC1, PC2, PC3, PC4, PC5',
    labSize = 3,
    returnPlot = FALSE)

pca_df <- as.data.frame(pca_result$rotated)

# Додаємо інформацію про condition та інші метадані
pca_df$condition <- colData_new$condition

# Використовуємо більше компонент для кластеризації, наприклад, перші 5
set.seed(1234)
km <- kmeans(x = pca_df[, 1:5], centers = 3, nstart = 25)  # Використовуємо 5 компонент і 25 запусків
# Додаємо результати кластеризації в датафрейм
pca_df$cluster <- as.factor(km$cluster)

# Візуалізація
pkmean <- autoplot(km, data = pca_df, frame = TRUE, frame.type = 'norm', size = 4, alpha = 0.6) + 
  aes(shape = condition, color = cluster) + 
  scale_shape_manual(values = c('0Gy' = 16, '2Gy' = 17, '4Gy' = 18)) + 
  labs(subtitle = 'with 5 PCs',
       color = 'cluster', shape = 'Condition') +
  coord_equal() + 
  theme_minimal(base_size = 18)+
  theme(legend.position = 'right', # Позиція легенди
        panel.border = element_rect(color = 'black', fill = NA, size = 1))  # Додання рамки


# Додаємо заголовки до кожного графіка
pscree <- pscree + ggtitle("A Scree plot") + theme(plot.title = element_text(face = "bold", hjust = -0.01, size = 18))
ploadings <- ploadings + ggtitle("B Loading plot") + theme(plot.title = element_text(face = "bold", hjust = -0.01, size = 18))
pbiplot <- pbiplot + ggtitle("C Bi-plot") + theme(plot.title = element_text(face = "bold", hjust = -0.01, size = 18))
pkmean <- pkmean + ggtitle("D K-means Clusters") + theme(plot.title = element_text(face = "bold", hjust = -0.01, size = 18))
ppairs <- ppairs + ggtitle("E Pairs plot") + theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 18))

# Комбінування графіків у сітку
combined_plot <- (pscree | ploadings) / (pbiplot | pkmean) / ppairs


ggsave("./figs/combined_plot.png", combined_plot, width = 20, height = 30)

knitr::include_graphics("./figs/combined_plot.png")

Функціональний аналіз

Після проведення аналізу диференційної експресії було визначено роль цих генів в клітинах. Для цього було проведено аналіз збагачення за допомогою Enrichr

Функціональний аналіз Enrichr

Узагальнене

knitr::include_graphics("./figs/enrichr_results_all_up.png")

2Gy 0Gy

knitr::include_graphics("./figs/enrichr_results_2Gy_0Gy_up.png")

4Gy 0Gy

knitr::include_graphics("./figs/enrichr_results_4Gy_0Gy_up.png")

4Gy 2Gy

knitr::include_graphics("./figs/enrichr_results_4Gy_2Gy_up.png")

Аналіз метаболічних шляхів за допомогою Pathview

Для того, щоб явно побачити роль генів з підвищеною та зі зниженою експресією у клітинах лінії А549 було проведено аналіз метаболічних шляхів циклу клітини, проліферації, апоптозу.

logFC <- path_4Gy_vs_0Gy$log2FoldChange
names(logFC) <- path_4Gy_vs_0Gy$entrezgene_id


pathview(gene.data = logFC, 
         pathway.id = "hsa05207", 
         species = "hsa", 
         #kegg.native = T
         ut.suffix = "Non-small cell lung cancer (специфічно для A549)")

pathview(gene.data = logFC, 
         pathway.id = "hsa04110", 
         species = "hsa", 
         #kegg.native = T
         ut.suffix = "Cell cycle")

pathview(gene.data = logFC, 
         pathway.id = "hsa04115", 
         species = "hsa", 
         #kegg.native = T
         ut.suffix = "p53 signaling pathway")

pathview(gene.data = logFC, 
         pathway.id = "hsa04725", 
         species = "hsa", 
         #kegg.native = T
         ut.suffix = "NAChR")

pathview(gene.data = logFC, 
         pathway.id = "hsa04210", 
         species = "hsa", 
         #kegg.native = T
         ut.suffix = "Apoptosis")

Канцерогенез

knitr::include_graphics("hsa05207.pathview.png")

Проліферація незважаючи на активацію RSK, загальний рівень клітинного поділу, ймовірно, знизиться через зменшення експресії генів E2F1 і CDC6.

Виживання клітин буде зменшене, оскільки зниження рівнів Survivin і BCL2 сприяє збільшенню апоптозу.

Апоптоз – ймовірно, підвищиться, через зниження активності антиапоптотичних білків (Survivin і BCL2).

Загальний результат: фінально ми очікуємо зниження проліферації і виживання клітин із збільшенням апоптозу, що може мати позитивний вплив у контексті боротьби з раковими клітинами.

Клітинний цикл

knitr::include_graphics("hsa04110.pathview.png")

Проліферація значно знизиться, через зниження експресії багатьох генів, що регулюють клітинний цикл, реплікацію ДНК і мітоз (CycA/B, CDK1/2, Plk1, CDC20, E2F1-3, MCM, Dbf4).

Виживання клітин частково збережеться, оскільки підвищена експресія Mdm2 може пригнічувати p53 і знижувати апоптоз.

Апоптоз ймовірно, знизиться через пригнічення p53 шляхом підвищення Mdm2, хоча інші механізми можуть впливати на загальну картину.

Загальний результат: значне зменшення клітинної проліферації через порушення ключових етапів клітинного циклу і реплікації, при цьому виживання клітин може залишитися відносно високим через інгібіцію апоптозу (Mdm2).

р53 сигнальний шлях

knitr::include_graphics("hsa04115.pathview.png")

Проліферація знизиться через активацію p21, зниження Cyclin B і B99, а також активність Sestrins, які інгібують проліферацію.

Виживання клітин зменшиться через зниження Bcl-2 і активацію зовнішнього апоптозу через Fas.

Апоптоз збільшиться через підвищену експресію Fas та PAG608, хоча MDM2 намагатиметься блокувати p53-опосередкований апоптоз.

Загальний висновок: Проліферація клітин значно знизиться через затримку в клітинному циклі та пригнічення mTOR сигналінгу (через Sestrins), а апоптоз буде активований через Fas та інші проапоптотичні фактори.

Апоптоз

knitr::include_graphics("hsa04210.pathview.png")

Під впливом X-rays відбувається значний запуск апоптозу через активацію зовнішнього шляху (FAS/FASL) та ослаблення антиапоптичних сигналів (Bcl-2, IAPXIAP, Colfeyrin). Це призводить до масової клітинної загибелі.

Холінергічний синапс

knitr::include_graphics("hsa04725.pathview.png")

Зниження рівня Bcl-2 робить клітини більш вразливими до апоптозу.

Але якщо аплікувати нікотином nAChR, то це може стимулювати активацію сигнальних шляхів PI3K/Akt та MAPK, що стимулюють антиапоптичні сигнали та виживання клітин.

Загальний висновок: Навіть при зниженій експресії Bcl-2, активація PI3K/Akt шляхів через нікотин може частково компенсувати цей дефіцит, стимулюючи клітинні механізми виживання. Однак, без належного рівня Bcl-2, захист клітини буде ослабленим.